/*
 * Projekt: GMU 2013 - Canny hranovy detektor OpenCL/CUDA
 * Autori: Pavel Nevrala - xnevra00 (xnevra00@stud.fit.vutbr.cz)
 *         Martin Jurak - xjurak02 (xjurak02@stud.fit.vutbr.cz)
 */

#pragma OPENCL EXTENSION cl_khr_byte_addressable_store : enable

#define PI 3.14159265359f

/*
 * Soubor obsahuje zdrojove kody pro kernely
 */

/*
 * Program provadi radkovou konvoluci masky se vstupnim obrazem. Okraje obrazu jsou reseny zrcadlenim obrazu kolem strany.
 * @input - vstupni obraz
 * @output - vystupni obraz
 * @width - sirka obrazu
 * @height - vyska obrazu
 * @mask - maska pro konvoluci
 * @maskSize - velikost masky
 */
__kernel void convolutionX ( __global float* input, __global float* output, uint width, uint height, __global float* mask, uint maskSize)
{
	int gx = (int) get_global_id(0); //defines column
	int gy = (int) get_global_id(1); //defines row

	int lx = get_local_id(0);
	int lw = get_local_size(0);
	int kw = 0;

	if( (gx < 0) || (gx > (width-1)) || (gy < 0) || (gy > (height-1) ))
		return;

	int actualPixel;
	int maskMid = maskSize/2; // stred masky-1

	// vysledny pixel po konvoluci
	float sum;
	float mask_weight;
	int offset;

	while(lx+kw < width){	

		//aktualni pixel 
		actualPixel = lx+kw + gy*width;
		sum = 0;

		// konvoluce s maskou pro dany pixel
		for(int i = 0; i<maskSize; i++){
			mask_weight = mask[i];
			// levy okraj obrazu		
			if((lx+kw - maskMid+i) < 0){
				offset = lx+kw - maskMid+i;
				sum += input[gy*width+(abs(offset)-1)] * mask_weight;
			}
			// pravy okraj obrazu
			else if((lx+kw + (i-maskMid)) > width-1){
				sum += input[actualPixel-(i+maskMid+1)] * mask_weight;
			}
			// zbytek obrazu
			else if((lx+kw < width) && ((lx+kw >= 0))){
				sum += input[actualPixel+i-maskMid] *mask_weight;
			}	
				
		}

		// vysledny pixel
		output[actualPixel]=sum;

		// pokud je sirka obrazu vetsi nez max pocet work-groupy, posune se na pixel = (aktualnipixel + sirka obrazu)
		kw += lw -1;		
	}
}

/*
 * Program provadi sloupcovou konvoluci masky se vstupnim obrazem. Okraje obrazu jsou reseny zrcadlenim obrazu kolem strany.
 * @input - vstupni obraz
 * @output - vystupni obraz
 * @width - sirka obrazu
 * @height - vyska obrazu
 * @mask - maska pro konvoluci
 * @maskSize - velikost masky
 */
__kernel void convolutionY ( __global float* input, __global float* output, uint width, uint height, __global float* mask, uint maskSize)
{
	int gx = (int) get_global_id(0);
	int gy = (int) get_global_id(1);

	int ly = get_local_id(1);
	int lh = get_local_size(1);
	int kh = 0;

	if( (gx < 0) || (gx > (width-1)) || (gy < 0) || (gy > (height-1) ))
		return;

	int actualPixel;
	int maskMid = maskSize/2; // stred masky-1

	// vysledny pixel po konvoluci
	float sum;
	float mask_weight;
	int offset;
	// konvoluce s maskou pro dany pixel
	while(ly+kh < height){
		
		//aktualni pixel 
		actualPixel = (ly+kh)*width + gx;
		sum = 0;
		
		// konvoluce s maskou pro dany pixel
		for(int i = 0; i<maskSize; i++){
			mask_weight = mask[i];
			
			// horni okraj obrazu	
			if((ly+kh - maskMid+i) < 0){
				offset = (ly+kh - maskMid+i);
				sum += input[(abs(offset)-1) *width + gx] *mask_weight;

			// dolni okraj obrazu
			}else if((ly+kh + (i-maskMid+1)) > height-1){
				sum += input[(ly+kh-(i+maskMid+1))*width + gx]* mask_weight;

			// zbytek obrazu
			}else if((ly+kh < height)&&(ly+kh >= 0)){
				sum += input[(ly+kh+i-maskMid)*width + gx]* mask_weight;
			}
		}

		// vysledny pixel
		output[actualPixel]=sum;
		
		// pokud je vyska obrazu vetsi nez max pocet work-groupy, posune se na pixel = (aktualnipixel + vyska obrazu)	
		kh += lh-1;
	}	
}


// =======================================================================================================================
// =======================================================================================================================

/*
 * Provede vypocet velikosti a smeru gradientu v obraze
 * Velikost gradientu vypocet Eklidovou metrikou sqrt(GX^2,GY^2)
 * Smer gradientu je dan funkci atan2(GY,GX). Pote je uhel zaokrouhlen do jednoho z 4 uhlu 0, 45, 90, 135.
 * @inputX - gradient ve smeru X
 * @inputY - gradient ve smeru Y
 * @output - vystupni obraz
 * @edgesDirection - vystupni obraz pro smery gradientu
 * @width - sirka obrazu
 * @height - vyska obrazu 
 */
__kernel void sobelResult ( __global float* inputX, __global float* inputY, __global float* output, __global float* edgesDirection, uint width, uint height)
{
	int gx = get_global_id(0);
	int gy = get_global_id(1);

	int lx = get_local_id(0);
	int lw = get_local_size(0);

	if( (gx < 0) || (gx > (width-1)) || (gy < 0) || (gy > (height-1) ))
		return;

	float GX, GY;
	float direction;
	
	int kw = 0;
	while(lx+kw < width)
	{
		// aktualni pixel
		int actualPixel = gy*width + lx + kw;

		// hodnota GX a GY
		GX = inputX[actualPixel];
		GY = inputY[actualPixel];
		
		// velikost gradientu
		output[actualPixel] = hypot(GX, GY);		

		// smer hrany
		direction = atan2(GY,GX)/(PI)*180;
		if(direction < 0){
			direction = 360 + direction;
		}	
		
		//zaoktrouhleni na uhel 0		
		if(((direction <= 22.5f) && (direction >= 0.0f)) || ((direction <= 360.0f) && (direction >= 337.5f)) || ((direction <= 202.5f) && (direction >= 157.5f))){
			direction = 0.0f;
		}
		//zaoktrouhleni na uhel 45
		else if (((direction <= 67.5f) && (direction > 22.5f)) || ((direction <= 247.5f) && (direction > 202.5f))){
			direction = 45.0f;
		}
		//zaoktrouhleni na uhel 90
		else if (((direction <= 112.5f) && (direction > 67.5f)) || ((direction <= 292.5f) && (direction > 247.5f))){
			direction = 90.0f;
		}
		//zaoktrouhleni na uhel 135
		else if (((direction < 157.5f) && (direction > 112.5f)) || ((direction < 337.5f) && (direction > 292.5f))){
			direction = 135.0f;		
		}

		edgesDirection[actualPixel] = direction;	
		kw += lw-1;		
	}
}

// =======================================================================================================================
// =======================================================================================================================

/*
 * Provede vypocet lokalniho maxima v obrazu (ztenceni hran)
 * Podle smeru daneho pixelu jsou vybrany dva pixeli z osmikololi (napr. pro uhel 0 - predchozi a nasledujici pixel obrazu v radku)
 * a porovna jejich velikost gradientu. Pokud ma jeden z okolnich vetsi gradient je aktualni pixel potlacen nastaven na 0 jinak 
 * zustane zachovan.
 * @input - vstupni obraz s velikosti gradientu 
 * @edgesDirection - smery gradientu
 * @output - vystupni obraz
 * @width - sirka obrazu
 * @height - vyska obrazu 
 */
__kernel void nonMaxSup(__global float* input, __global float* edgesDirection,__global float* output, uint width, uint height)
{
	int gx = get_global_id(0);
	int gy = get_global_id(1);

	int lx = get_local_id(0);
	int lw = get_local_size(0);
	
	if( (gx < 0) || (gx > (width-1)) || (gy < 0) || (gy > (height-1) ))
		return;
	
	int kw = 0;
	float direction, center, first, second;
	int actualPixel;
	while(lx + kw < width)
	{	
		// aktualni pixel
		actualPixel = gy*width + lx + kw;
		// smer gradientu
		direction = edgesDirection[actualPixel];
		// hodnota gradientu
		center = input[actualPixel];		
		
		// dle smeru gradientu porovna pixely
		if(direction == 0 ){				
			if (lx+kw == 0){
				first = 0.0f;
				second = input[actualPixel+1]; 
			}else if (lx+kw == width-1){
				first = input[actualPixel-1]; 
				second = 0.0f;
			}else if ((lx+kw > 0) && (lx+kw < width-1)){
				first = input[actualPixel-1]; 
				second = input[actualPixel+1];
			}		
		}		
		else if(direction == 90){
			if (gy == 0){
				first = 0.0f;
				second = input[actualPixel+width]; 
			}else if (gy == height-1){
				first = input[actualPixel-width];
				second = 0.0f;
			}else if ((gy > 0) && (gy < height-1)){
				first = input[actualPixel-width];
				second = input[actualPixel+width];
			}
		
		}
		else if(direction == 45){
			if (((lx + kw == 0) && (gy == height-1)) || ((lx + kw == width-1) && (gy == 0))){
				first = 0.0f;
				second = 0.0f;
			}else if (((lx + kw < width) && (lx + kw > 0) && (gy == height-1)) || ((lx + kw == width-1 ) && (gy < height) && (gy > 0))){
				first = input[actualPixel-width-1];
				second = 0.0f;
			}else if (((lx + kw < width) && (gy == 0)) || ((lx + kw == 0) && (gy < height))){
				first = 0.0f;
				second = input[actualPixel+width+1];
			}else if ((lx + kw > 0) && (lx + kw < width) && (gy > 0) && (gy < height-1)){
				first = input[actualPixel-width-1];
				second = input[actualPixel+width+1];
			}
		}else if(direction == 135){
			if (((lx + kw == 0) && (gy == 0)) || ((lx + kw == width-1) && (gy == height-1))){
				first = 0.0f;
				second = 0.0f;
			}else if (((lx + kw == width-1) && (gy < height-1) && (gy > 0)) || ((lx + kw > 0) && (lx + kw < width) && (gy == 0))){
				first = input[actualPixel+ width-1];
				second = 0.0f;
			}else if (((lx + kw == 0) && (gy < height) && (gy > 0)) || ((lx + kw >= 0) && (lx + kw < width-1) && (gy == height-1))){
				first = 0.0f;
				second = input[actualPixel - width+1];
			}else if ((lx + kw > 0) && (lx + kw < width-1) && (gy > 0) && (gy < height-1)){
				first = input[actualPixel+ width-1];
				second = input[actualPixel- width+1];
			}
		}
		
		// porovnani
		if((center < first) || (center < second) || (((lx+kw== width-1) && (lx+kw == 0)) && ((gy== height-1) && (gy == 0))) ){
			output[actualPixel] = 0.0f;
		}			
		else{
			output[actualPixel] = center;
		}

		kw += lw - 1;
	}

}

// =======================================================================================================================
// =======================================================================================================================

/*
 * Provede doublethresholding s hysterezi
 * Prochazi pixel po pixelu. Pokud se jeden o pixel na silne hrane, zjisti smer gradientu v tomto bode
 * a zkouma pixel v tomto smeru. Pokud se jedna o slabou hranu (mezi dolnim a hornim prahem) oznaciji ji jako silnou.
 * Ostatni pixeli, ktere nesjou napojene na silne hrany jsou nastaveny na intezitu 0.
 * @input - vstupni obraz s velikosti gradientu 
 * @edgesDirection - smery gradientu
 * @output - vystupni obraz
 * @thrLow - vystupni obraz
 * @thrHigh - vystupni obraz
 * @width - sirka obrazu
 * @height - vyska obrazu 
 */
__kernel void doubleThresholding( __global float* input,  __global float* edgesDirection, __global float* output, uchar thrLow, uchar thrHigh, uint width, uint height)
{
	int gx = get_global_id(0);
	int gy = get_global_id(1);		
	
	int lx = get_local_id(0);
	int lw = get_local_size(0);

	if( (gx < 0) || (gx > (width-1)) || (gy < 0) || (gy > (height-1) ))
		return;

	float pixel;
	int actualPixel = gx + gy*width;
	pixel = input[actualPixel];
	float direction = edgesDirection[actualPixel];
	
	int new_x = gx;
	int new_y = gy;
	if(pixel >= convert_float(thrHigh)){

		// pixel lezi na silne hrane
		// podle smaru gradientu prochazi okolni pixely
		if(direction == 0 ){				
			while((new_x >= 0) && (new_x < width) && (edgesDirection[new_x + new_y*width] == 0) && (input[new_x + new_y*width]>thrLow)){
				output[new_x + new_y*width]=255;
				new_x = new_x+1;
                new_y = new_y;					                   
			}
		}
		else if(direction == 90){
			while((new_y >= 0) && (new_y < height) && (edgesDirection[new_x + new_y*width] == 90) && (input[new_x + new_y*width]>thrLow)){
				output[new_x + new_y*width]=255;	
				new_x = new_x;
                new_y = new_y+1;					                  
			}				
		}
		else if(direction == 45){
			while((new_y >= 0) && (new_y < height) && (new_x >= 0) && (new_x < width) && (edgesDirection[new_x + new_y*width] == 45) && (input[new_x + new_y*width]>thrLow)){
				output[new_x + new_y*width]=255;
				new_x = new_x+1;
                new_y = new_y-1;					               
			}
		}else if(direction == 135){
			while((new_y >= 0) && (new_y < height) && (new_x >= 0) && (new_x < width) && (edgesDirection[new_x + new_y*width] == 135) && (input[new_x + new_y*width]>thrLow)){
				output[new_x + new_y*width]=255;
				new_x = new_x+1;
                new_y = new_y+1;                  
			}
		}			
	}
	else
		output[actualPixel] = 0;
}